% code to compute MT trfn from raw ts files 
% Support RTS files created by GMS05 
% usage -> mtprocn 
% uses a suit of sub-functions.

% references-----------
% 1) Sims, 1971 (ls impedance) 2) Dekker, (predicted coherency) 1984 
% 3) Gamble (rr) , 1979 (rr)	4)ProcMT manual (rts file structue, biv error)
% 5) Manoj & Nagarajan, 2003 (ANN selection)
% 6) Chave et al 1989 (Jack Knife) 7) Ritter et al 1998 (Robust processing)
% developement history-----
% 8/2000 - compute cal. table according to gain settings etc;
% 11.7.01 -computes cross and auto spectra & average it 
% 12.7 - MT app rest & phse are now computed
% 8/02 - e-pred coherency, and Biv error are calculated
% 10.2002 - rr processing introduced
% 12.2002 - rr error is being introduced
% 04.03 -- Jack knife impedance & variance
% 05.03.03 ann selection for ts segments
% 10.03.03 robspm Robust averaging for SPMs
% 15.03.03 normalize SPMs wrt power in horizontal magnetic fileds
% Robust band average and 10.6.3
% Robust section average using Jack knife 10.6.3


% ----known bugs & fixings
% phase distortion for 2nd and 3rd band- 05.07.2002
% OK for band1 & 4 on 07.07.2002
% 10.02 fixed phase distortion (problem was with complex transpose)
% 11.02 band2 & band1 is recorded as slices of 1024 length
%       this is not taken care of (nfft & overalp) 
% 11.02 computed error for rr is still local one
% 12.02 band1 phase distortion - problem with bd1 sys trfn ?

% ------future projects---
% ann despiker			depiker attempted bu with partial success 31.5.3
% phase rotation for time synchronization problems (?)  % Jones 1983 gave an equation for this but not yet implimented
% write edi to be completed			% not yet
% ats plugg-in -> for ADU06 ts		% not yet implimented
% latest date 02.07.2003

[fid,er,file,path1] = getbfile('c:\manoj\ngri-dell01\*.rts');

if fid == -1 | er < 1,
 	return;
end;

disp('Rts-in -> reading time series');
[TS,fid] = rts_in(fid); % reading time series file (rts format - metronix)
fclose(fid);

%%%%%%%------------------------------ be careful delete te following LIN
TS.head.fs = 1;


ProcDef = SetDefault(TS.head.band); % setting defaults. ProcDef is a structured array
ProcDef.file1=file;

if ProcDef.ref == 1, % rr processing
   disp('Show me the reference time series !');
   [fid,er,file] = getbfile('c:\manoj\ngri-dell01\*.rts');
   if fid == -1 | er < 1,
 	return;
end;

disp('Rts-in -> reading reference time series');
  	[TSr,fid] = rts_in(fid);
  	fclose(fid);
  	disp('RrCorrectTS -> Checking overlap, trimming offshoots.. ');
  	[TS,er,len] = RrCorrectTS(TS,TSr,ProcDef); % check for overlap, trim extra lengths, pack the reference field to the local filed
  
  	if er<1,
	   fprintf('Error file mismatch\n');
	   er=-1;
   return;
   else,
   	ProcDef.file2=file;
   	clear TSr;
   end;
end;

ProcDef.len = floor(length(TS.matrix)/ProcDef.block);

ProcDef.SelStacks = ones([1 ProcDef.len]);  

fprintf('File %s, Band %d, Stacks %d Nfft %d Overlap %d\n', ProcDef.file1,ProcDef.nband,ProcDef.len,ProcDef.block,...
   ProcDef.overlap);
fprintf('Processing parameters: Section Averaging %s, Band Averaging %s\n',ProcDef.Method,ProcDef.BandAv);

if ProcDef.nband == 4 & ProcDef.DoANN == 1,
	ProcDef = SelANN(TS,ProcDef);
end;                                               
                                               

if ProcDef.BandAv == 'R', %robust band averaging 
   [SPMatrix] = GetRobSpectra(TS,ProcDef);
else, % Parzen band averagibg
   [SPMatrix] = GetSpectra(TS,ProcDef);
end;
   

temp = size(SPMatrix); % This is required, as number of overlapped ts sections may not be just 2*len-1
ProcDef.len = temp(2);

ProcDef.Weights = ones([ProcDef.nfrq ProcDef.len]); % Weights are modified in robust averaging routine(s)

SPMatrix = ScaleSPM(SPMatrix,ProcDef);

%disp('NormSPM-> Normalize the SPMs wrt power in horz. mag fields');
%SPMatrix = NormSPM(SPMatrix,ProcDef); % Normalize the SPMs wrt power in horz. mag fields


switch ProcDef.Method,
case 'LS', % Least Square
   disp('AverageSPM (Non robust)-> Averaging  spectra sets....');
   SPM2  = AverageSPM(SPMatrix,ProcDef);
	disp('TransFn   - > Computing impedances....');
   tf2=TransFn(SPM2,ProcDef);
case 'RB', % Robust M
   [tf2,ProcDef,SPM2] = robspm(SPMatrix,ProcDef);
case 'CT', % Coherency Thrashold
   [SPM2,ProcDef]= cohspm(SPMatrix,ProcDef,2);
   [tf2]=tf_dn(SPM2);
   disp('Warning Errors are not computed\n');
case 'CO', % Preferential Coherency Stacking
	[SPM2,ProcDef]= cohspm(SPMatrix,ProcDef,1);
   [tf2]=tf_dn(SPM2);
   disp('Warning Errors are not computed\n');
case 'JK', % Jack knife impedance & error 
   [tf2,dof]=jknf(SPMatrix,ProcDef);
end;

   

figure(1);
disp('DispResPh - > Dispalying App Res,Ph & Errors..');

[Rhxy Rhyx Phxy Phyx] = dispresph(tf2,ProcDef);
dispresph(tf2,ProcDef);
if ProcDef.ref == 1,
   tf1 = tf_r(SPM2);
      figure(2);
      [Rhxy Rhyx Phxy Phyx] = dispresph(tf1,ProcDef);
      disp('Warning remote reference variance computed is not correct');
      disp('Local variance is shown');
end;


if ProcDef.WriteEdi == 1,
   fprintf('Edi-out -> writing out edi\n');
   ll=findstr(file,'.');
   file1=file(1:ll-1);
   file1(ll:ll+3)='.edi';
   fid = fopen([path1 file1],'wt');
   if fid == -1,
      fprintf('Unable to open %s error in edi_out\n',[path file1]);
      return;
   end;
   [fid,er]=edi_out(SPM2,fid);
   fclose(fid);
end;
